install.packages(c("devtools", "knitr", "rmarkdown"))
library(devtools)
install_github(c("MRCIEU/TwoSampleMR","MRCIEU/MRInstruments"))
install_github("WSpiller/MRPracticals",build_opts = c("--no-resave-data", "--no-manual"),build_vignettes = TRUE)

vignette("MRBase")
vignette("RadialMR")
vignette("MVMR Tutorial")

#心脏骤停MR，21种
setwd("D:\\ProgramFiles\\R\\Work\\CAMR")
library(TwoSampleMR)

load("ca.rdata")

exp_list <- read.delim("PMID37106081.tsv", header = F)
colnames(exp_list) <- c("ID", "exp","exposure")
head(exp_list)

exps <- lapply(exp_list$ID, function (x){
  aa<-read.delim(paste0("./exp/", x, "_buildGRCh37.tsv"))
  ab <- aa[aa$p_value < 5e-08,]
  })
names(exps) <- exp_list$exp
rows <- lapply(exps, function (x) {nrow(x)})
exps <- exps[rows != 0]


#用广州写的这个函数添加eaf。https://github.com/HaobinZhou/Get_MR
library(data.table)
get_eaf_from_1000G<-function(dat,path,type="exposure"){
  
  corrected_eaf_expo<-function(data_MAF){
    effect=data_MAF$effect_allele.exposure
    other=data_MAF$other_allele.exposure
    A1=data_MAF$A1
    A2=data_MAF$A2
    MAF_num=data_MAF$MAF
    EAF_num=1-MAF_num
    
    harna<-is.na(data_MAF$A1)
    harna<-data_MAF$SNP[which(harna==T)]
    
    cor1<-which(data_MAF$effect_allele.exposure !=data_MAF$A1)
    
    
    data_MAF$eaf.exposure=data_MAF$MAF
    data_MAF$type="raw"
    data_MAF$eaf.exposure[cor1]=EAF_num[cor1]
    data_MAF$type[cor1]="corrected"
    cor2<-which(data_MAF$other_allele.exposure ==data_MAF$A1)
    cor21<-setdiff(cor2,cor1)
    cor12<-setdiff(cor1,cor2)
    error<-c(cor12,cor21)
    data_MAF$eaf.exposure[error]=NA
    data_MAF$type[error]="error"
    
    data_MAF<-list(data_MAF=data_MAF,cor1=cor1,harna=harna,error=error)
    
    return(data_MAF)
    
  }
  
  corrected_eaf_out<-function(data_MAF){
    effect=data_MAF$effect_allele.outcome
    other=data_MAF$other_allele.outcome
    A1=data_MAF$A1
    A2=data_MAF$A2
    MAF_num=data_MAF$MAF
    EAF_num=1-MAF_num
    
    harna<-is.na(data_MAF$A1)
    harna<-data_MAF$SNP[which(harna==T)]
    
    cor1<-which(data_MAF$effect_allele.outcome !=data_MAF$A1)
    
    
    data_MAF$eaf.outcome=data_MAF$MAF
    data_MAF$type="raw"
    data_MAF$eaf.outcome[cor1]=EAF_num[cor1]
    data_MAF$type[cor1]="corrected"
    cor2<-which(data_MAF$other_allele.outcome ==data_MAF$A1)
    cor21<-setdiff(cor2,cor1)
    cor12<-setdiff(cor1,cor2)
    error<-c(cor12,cor21)
    data_MAF$eaf.outcome[error]=NA
    data_MAF$type[error]="error"
    
    data_MAF<-list(data_MAF=data_MAF,cor1=cor1,harna=harna,error=error)
    
    return(data_MAF)
    
  }
  
  if(type=="exposure" && (is.na(dat$eaf.exposure[1])==T || is.null(dat$eaf.exposure)==T)){
    r<-nrow(dat)
    
    setwd(path)
    MAF<-fread("fileFrequency.frq",header = T)
    
    dat<-merge(dat,MAF,by.x = "SNP",by.y = "SNP",all.x = T)
    
    dat<-corrected_eaf_expo(dat)
    
    cor1<-dat$cor1
    
    harna<-dat$harna
    
    error<-dat$error
    
    dat<-dat$data_MAF
    
    print(paste0("一共有",(r-length(harna)-length(error)),"个SNP成功匹配EAF,占比",(r-length(harna)-length(error))/r*100,"%"))
    
    print(paste0("一共有",length(cor1),"个SNP是major allele，EAF被计算为1-MAF,在成功匹配数目中占比",length(cor1)/(r-length(harna)-length(error))*100,"%"))
    
    print(paste0("一共有",length(harna),"个SNP在1000G中找不到，占比",length(harna)/r*100,"%"))
    
    print(paste0("一共有",length(error),"个SNP在输入数据与1000G中效应列与参照列，将剔除eaf，占比",length(error)/r*100,"%"))
    
    print("输出数据中的type列说明：")
    print("raw：EAF直接等于1000G里的MAF数值，因为效应列是minor allele")
    print('corrected：EAF等于1000G中1-MAF，因为效应列是major allele')
    print("error：输入数据与1000G里面提供的数据完全不一致，比如这个SNP输入的效应列是C，参照列是G，但是1000G提供的是A-T，这种情况下，EAF会被清空（NA），当成匹配失败")
    
    return(dat)
  }
  
  if(type=="outcome" && (is.na(dat$eaf.outcome[1])==T || is.null(dat$eaf.outcome)==T)){
    r<-nrow(dat)
    
    setwd(path)
    MAF<-fread("fileFrequency.frq",header = T)
    
    dat<-merge(dat,MAF,by.x = "SNP",by.y = "SNP",all.x = T)
    
    dat<-corrected_eaf_out(dat)
    
    cor1<-dat$cor1
    
    harna<-dat$harna
    
    error<-dat$error
    
    dat<-dat$data_MAF
    
    print(paste0("一共有",(r-length(harna)-length(error)),"个SNP成功匹配EAF,占比",(r-length(harna)-length(error))/r*100,"%"))
    
    print(paste0("一共有",length(cor1),"个SNP是major allele，EAF被计算为1-MAF,在成功匹配数目中占比",length(cor1)/(r-length(harna)-length(error))*100,"%"))
    
    print(paste0("一共有",length(harna),"个SNP在1000G找不到，占比",length(harna)/r*100,"%"))
    
    print(paste0("一共有",length(error),"个SNP在输入数据与1000G中效应列与参照列，将剔除eaf，占比",length(error)/r*100,"%"))
    
    print("输出数据中的type列说明：")
    print("raw：EAF直接等于1000G里的MAF数值，因为效应列是minor allele")
    print('corrected：EAF等于1000G中1-MAF，因为效应列是major allele')
    print("error：输入数据与1000G里面提供的数据完全不一致，比如这个SNP输入的效应列是C，参照列是G，但是1000G提供的是A-T，这种情况下，EAF会被清空（NA），当成匹配失败")
    
    return(dat)
  }
  else{return(dat)}
}
exps_cl <- lapply(names(exps), function (x) {
  exposure_data <- format_data(exps[[x]],snp_col = "variant_id", beta_col = "beta",
                             se_col = "standard_error", effect_allele_col = "effect_allele",
                             other_allele_col = "other_allele", pval_col = "p_value",
                             ncase_col = "N")
  exposure_data <- get_eaf_from_1000G(exposure_data, ".\\", type = "exposure")
  exposure_data <- clump_data(exposure_data, clump_r2 = 0.01, pop = "EUR")
  exposure_data$exposure <-x
  exposure_data$id.exposure <- x
  return(exposure_data)
})

exp_mix <- do.call(rbind, exps_cl)
out_data <- read_outcome_data(snps = exp_mix$SNP, filename = "finngen_R10_I9_CARDARR.gz",
                               sep = "\t", snp_col = "rsids",beta_col = "beta",
                               se_col = "sebeta", effect_allele_col = "alt", other_allele_col = "ref",
                               pval_col = "pval", gene_col = "nearest_genes", chr_col = "X.chrom",
                               pos_col = "pos")
out <- get_eaf_from_1000G(out_data, ".\\", type = "outcome")
data <- harmonise_data(exp_mix, out)
res <- mr(data, method_list = c("mr_wald_ratio", "mr_ivw", "mr_egger_regression",
                                "mr_weighted_median", "mr_weighted_mode"))
he <- mr_heterogeneity(data)
ho <- mr_pleiotropy_test(data)
p1 <- mr_scatter_plot(res, data)
write.table(he, file = "he.txt", sep = ",", row.names = F, col.names = T, quote = F )
write.table(ho, file = "ho.txt", sep = ",", row.names = F, col.names = T, quote = F )
pdf("ldl1.pdf", height = 8, width = 10)
print(p1[11])
dev.off()
pdf("sbp1.pdf", height = 8, width = 10)
print(p1[18])
dev.off()
res_single <- mr_singlesnp(data, all_method = c("mr_ivw", "mr_egger_regression",
                                                "mr_weighted_median", "mr_weighted_mode"))
p2 <- mr_forest_plot(res_single)
p3 <- mr_funnel_plot(res_single)

res_loo <- mr_leaveoneout(data)
p4 <- mr_leaveoneout_plot(res_loo)

pdf("ldl2.pdf", height = 18, width = 10)
print(p2[11])
dev.off()
pdf("sbp2.pdf", height = 16, width = 10)
print(p2[18])
dev.off()
pdf("ldl3.pdf", height = 8, width = 10)
print(p3[11])
dev.off()
pdf("sbp3.pdf", height = 8, width = 10)
print(p3[18])
dev.off()
pdf("ldl4.pdf", height = 18, width = 10)
print(p4[11])
dev.off()
pdf("sbp4.pdf", height = 16, width = 10)
print(p4[18])
dev.off()

res$pas <- paste0(res$exposure, " ", res$method)
res$pas <- factor(res$pas, levels = rev(res$pas))
res$or <- exp(res$b)
res$ci_u <- exp(res$b + 1.96*res$se)
res$ci_d <- exp(res$b - 1.96*res$se)

write.table(res, "res.csv", row.names = F, col.names = T, sep = ",", quote = F)

p_sl <- ggplot(res, aes(x = or, y = pas))+
  geom_rect(aes(xmin = -Inf, xmax = Inf, ymin = which(pas == pas) - 0.25,
                     ymax = which(pas == pas) + 0.25), fill = "grey95")+
  geom_point(aes(col = exposure))+
  scale_y_discrete(labels = c("Wald ratio",rep(c("Weighted mode","Weighted median",
                                                 "MR Egger","IVW"),times = 19)))+
  geom_errorbarh(aes(xmax = ifelse(ci_u >3, 2.97 , ci_u), 
                     xmin = ci_d, col = exposure), height = 0.5)+
  geom_segment(aes(x = ifelse(ci_u >3, 2.97, 5), xend = ifelse(ci_u >3, 3, 6), y = pas, yend = pas, colour = exposure),
               arrow = arrow(type = "closed", length = unit(1.5, "mm")))+
  geom_vline(xintercept = 1)+
  scale_x_continuous(limits = c(0, 3))+
  theme_classic()+
  theme(legend.position = "none",
        axis.ticks = element_blank(),
        axis.line.y = element_blank(),
        axis.text.y = element_text(color = "black"))+
  labs(y = "", x = "")+
  geom_text(aes(x = 3, y = pas, label = formatC(pval, format = "e", digits = 2)),
                size = 3, color = ifelse(res$pval < 0.002, "red", "black"))
  
res$exposure <- factor(res$exposure, levels = rev(unique(res$exposure)))
p_na <- ggplot(res, aes(x = outcome, y = exposure, label = exposure))+
  geom_text(size = 3, color = "black")+
  theme(panel.grid.major = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white"))+
  labs(x = "Exposures")+
  scale_x_discrete(position = "top")

res$nsnp <- as.character(res$nsnp)
p_snp <- ggplot(res, aes(x = outcome, y = exposure, fill = nsnp))+
  geom_text(aes(label = nsnp),size = 3, color = "black")+
  theme(panel.grid.major = element_blank(),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y = element_blank(),
        panel.background = element_rect(fill = "white"))+
  labs(x = "Number of SNPs")+
  scale_x_discrete(position = "top")  

pdf("mr.s.pdf", height = 10, width = 10)
p_na + p_sl
dev.off()
